Summary. — Granular media such as sand and sugar are ubiquitous in nature and 
industry but are less well understood than fluids or solids. We consider the behavior 
of rapid granular flows where the transfer of momenta by collisions dominates. The 
physics is quite different for the opposite limit of static or slowly moving grains (e.g., 
sand piles). To gain understanding of granular flows we consider two problems that 
have been investigated with experiments, particle simulations and hydrodynamic 
theory: vertically oscillating granular layers and flow past an obstacle. Oscillating 
granular layers spontaneously form spatial patterns when the container acceleration 
amplitude exceeds a critical value, about 2.5 times the gravitational acceleration. 
Simulations with hard spheres that conserve linear momentum and dissipate energy 
in collisions are in qualitative accord with some but not all aspects of the observed 
patterns. It is necessary to include friction and angular momentum conservation in 
the simulations to achieve quantitative accord with observations. 
The applicability of a hydrodynamic theory to granular flows is not obvious because 
for typical conditions the particle mean free path is comparable to the length scale 
over which velocity and density fields change; hence there is not the separation of 
scales needed to justify a hydrodynamic approach. However, we show that Navier- 
Stokes-like equations describe well the density and temperature fields in vertically 
oscillating layers, even though this system is far from isothermal, incompressible 
fluid dynamics. 

The second problem examined, flow past an obstacle, is also described well by par- 
ticle simulations for frictional dissipative particles, but continuum simulations are 
only in qualitative accord with the observations. We show that continuum theory 
fails because it does not include friction between particles. Finally, we discuss how 
shock waves are common in granular flows since the the speed of sound (pressure 
waves) in a granular gas is typically only a few centimeters/second, while mean 
flow speeds are typically meters/second. Comparison of shocks in granular experi- 
ments and simulations is made with shocks in ordinary gases. Much more research 
is needed to understand how shocks evolve in granular flows. 
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1. Introduction 

Granular materials include sand, sugar, crushed coal, cereals, pills, cosmetics, and 
asteroids. The transport, mixing, and segregation of granular materials is important in 
the pharmaceutical, mining, agricultural, metal, food, and energy industries. Even in 
the chemical industry, the majority of the products are in granular rather than liquid 
form [1]. Thus a large engineering literature has developed on "powders and partic- 
ulates". However, a basic understanding of the physical mechanisms underlying the 
collective behavior of particles in a granular medium is lacking. Granular media can 
exhibit both solid and fluid properties (e.g., one can walk on a beach or pour the sand 
from a bucket), but granular media are less well understood than solids and fluids. While 
fluids are processed in industry with high efficiency, the efficiency of handling (crushing, 
mixing, separation) of granular materials is estimated to be well below optimum [1]. 

The scientific study of granular systems has a long history, including a discussion by 
Galileo in his Dialogues and an 1831 study by Faraday of convective motion of grains 
in heaps in vertically oscillated granular layers [2]. In recent years there has been a 
resurgence of interest in granular media among physicists, thanks to Pierre-Giles de 
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Genncs, who recommended in the early 1980s to young French scientists that granular 
matter was an interesting subject worthy of study [3-6]. In the past two decades there 
has been an explosion of interest in granular systems. A search in INSPEC on the word 
"granular" followed by "system, medium, matter, flow, or gas" yields 20 papers in the 
three year period 1980-82, while a decade later the number of papers in a three year period 
jumped to 112, and in 2000-2002 there were 691 papers on the subject. The growth is not 
so dramatic for literature searches using terms often used in industry (powder, particles, 
particulates), but most of that work is empirical, and much of that literature concerns 
pastes, soil, and fine particles (powders). We make no attempt to review the enormous 
literature on granular/particulate matter but list a few reviews [7-13], books [14-18], and 
several journals: Powder Technology, Granular Matter, Particle Science and Technology, 
and Advanced Powder Technology. 

The key property distinguishing collisions of granular particles from collisions of atoms 
in an ordinary gas is dissipation: collisions between macroscopic grains are inelastic - 
in the absence of external forcing, the particles in a granular medium all come to rest. 
A measure of this dissipation is the coefficient of restitution e, which is the ratio of the 
relative normal velocity of two particles after a collision to the relative normal velocity 
before collision. Some representative values of e for a relative normal velocity of 0.1 m/s 
are 0.96 for hardened bronze, 0.85 for aluminum, and 0.3 for lead [19]; the values depend 
on the particle bulk and surface properties in a complicated way (see, e.g., [20]). 

This chapter concerns rapid granular flows, which is called the collisional regime to 
distinguish it from the quasi-static regime where particles are at rest or nearly so. Much 
of the granular literature concerns the quasi-static regime where inertia is not significant 
and chains of particles in contact bear most of the load. Understanding the development 
and evolution of these force chains and the role of steric hindrance is often the focus 
of the research on sand piles and other granular systems at rest or nearly at rest. The 
collisional and quasi-static regimes are each difficult, but the intermediate regime with 
some particles moving rapidly while others are at rest is even more difficult. For an 
example, see the contribution in this book on the formation of craters in a granular 
medium [21,22]. 

We will consider situations where only contact forces are important. For small parti- 
cles (less than«50 /um), electrostatic and van der Waals forces become important. Fur- 
ther, air friction can be significant, but for particles greater than about 1 mm in diameter, 
air friction is often negligible if the velocity is not too large. 

An oft-studied example of a granular medium in the collisional regime is a collection of 
particles in a vertically oscillating container. Section 2 describes spatial patterns formed 
by such a system. Since the number of particles in an experiment can be small (less 
than one million) , Newton's laws for the motion of these particles can be directly imple- 
mented on a Personal Computer. Section 3 discusses such Molecular Dynamics (MD) 
simulations for hard, spherical particles that are characterized by a restitution coefficient 
and a frictional coefficient. Most theoretical analyses of granular flows examine fric- 
tionless (smooth) inelastic spheres, but there exist no frictionless macroscopic particles, 
just as there are no elastic particles. Molecular dynamics simulations show that realistic 
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models of rapid granular motions must include friction - friction provides another mode 
of dissipation and also results in reduced grain mobility and a higher particle collision 
rate [23]. Molecular dynamics simulations including friction describe experimental ob- 
servations on oscillating granular media very well, as we shall describe in Section 3, but 
simulations without friction fail to capture even qualitatively some important aspects of 
the observations. 

Studies of oscillated granular layers have revealed localized structures, "oscillons", 
which are stable for a range of container oscillation frequencies and amplitudes. Section 
4 describes the properties of oscillons and shows how they can be considered as the basic 
building blocks of some extended spatial patterns. 

The continuum approach to granular flows is introduced in Section 5, where equations 
of motion derived by Jenkins and Richman for a dilute dissipative fluid are presented. 
The Jenkins-Richman equations are similar to the Navier-Stokes equations for a fluid 
but are modified to include dissipative effects arising from collisions between the grains. 
However, the equations do not include the effect of friction in grain-grain collisions or 
grain-container collisions. As will be discussed, the application of continuum equations 
to granular flows must be considered with caution because, unlike an ordinary fluid, in 
granular flows there is not a large separation between the microscopic and macroscopic 
length and time scales. 

In a granular fluid the speed of propagation of pressure fluctuations (the sound speed) 
is only a few centimeters/sec, which is very small compared to the 330 m/s sound speed in 
air. (The granular sound speed here refers to particles in vacuum; there is no air and the 
only particles are the dissipative grains.) Hence the average streaming speed of granular 
flows generally exceeds the sound speed in the flow, and shock waves form whenever 
the flow encounters an obstacle. These shock waves are the subject of Section 6, where 
observations from experiments are compared with molecular dynamics and continuum 
simulations. 

We will conclude with a discussion in Section 7 that mentions some of the many issues 
open in the understanding of systems of dissipative particles. 

2. Patterns in vertically oscillated granular layers 

Consider a layer of non-cohesive grains in a container oscillating sinusoidally with 
frequency / and dimensionless acceleration amplitude T = 4ir 2 f 2 A/ g, where 2A is the 
peak-to-peak amplitude of the displacement of the container and g is the gravitational 
acceleration. We consider patterns that arise spontaneously, not from sidewall forcing 
or from interstitial gas, but from correlations induced by multiple collisions between the 
grains and between the grains and the container bottom. To minimize sidewall effects, 
the horizontal dimensions of the container are made large compared to the layer depth 
h. The layer is illuminated from the side, as fig. ^ illustrates [24,25]. 

For r > 1, on each cycle the layer loses contact with the container, flies in the air, 
and then collides with the container. However, the layer remains compact and flat until 
r w 2.5, where a standing wave pattern spontaneously forms. A square pattern forms at 
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Fig. 1. - (a) Electromechanical shaker system for studying pattern formation in vertically oscil- 
lating layers, (b) Spatial patterns are illuminated from the side by light incident at low angles. 
When viewed from above, high regions are bright and low regions are dark. From [24]. 



low frequencies and a stripe pattern at high frequencies, as figs. Eta) and (b) respectively 
illustrate [26,27]. The pattern is subharmonic, repeating every 2r, where r = 1// is 
the container oscillation period; thus a ridge in a striped pattern at an instant of time 
becomes a valley one container oscillation period later. 

A heuristic argument for the critical value of T for the onset of instability of a flat 
oscillating layer was given in [25] . The authors argued that a flat layer becomes unstable 
when the collision occurs at the plate's lowest point. Then a time r/2 is taken for the 
layer to free fall from its highest point through a distance 2A to collide with the plate, so 
that 2A — i(ji(r/2) 2 = \gr 2 , which gives T c = 2.5, in accord with the observed critical 
acceleration for the onset of patterns. 

The onset of patterns (squares or stripes) at T w 2.5 is quite robust, independent of 
layer depth (depths from about a monolayer to about 25particlediametersa), container 
shape, and particle properties (size, restitution coefficient, surface roughness, material). 
Most studies of patterns have been conducted for particles of uniform size, but studies 
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Fig. 2. - Patterns in oscillating granular layers: (a) stripes, (b) squares, and (c) hexagons (with 
different phases on the left and right). The patterns oscillate at //2, where / is the container 
oscillation frequency. Here F, /, and h/a (ratio of the depth of the layer at rest to the diameter 
of the particles) are given by (a) 3.3, 67 Hz, 7, (b) 2.9, 25 Hz, 4, and (c) and (d) 4.0, 67 Hz, 
7. In (a) and (c) the diameter of the container is 770a and in (b) the container is square, 
1100(7 x llOOcr. The particles are bronze spheres 0.165 mm diameter, (a) and (c) are from [26] 
and (b) is from [27]. 



with a range of particle sizes have found that the onset remains sharp for size distributions 
ranging up to about 30% [28]. Square and stripe patterns have also been observed to 
form at T 2.5 for irregular particles such as rice grains and grass seed [28]. 

The phase diagram (fig. [3J shows the stability regions for different patterns as a 
function of T and dimensionless frequency /*, where /* = \f~hfg and h is the depth of 
the layer at rest [29]. The transitions are well defined and are only weakly dependent 
on /*. Except for the transition from a flat layer to squares, the hysteresis is small or 
perhaps zero. We know of no argument for the square to stripe transition as a function 
of frequency with T fixed, but we note that the transition has been found to occur at 
/* « 1/3 [24]. 

Some insight into the dynamics can be gained from consideration of the one-dimensional 
motion of a single completely inelastic ball on a vertically oscillating plate [30] . Such a 
model cannot of course describe the 2D spatial patterns that form in granular layers on 
an oscillating plate, but it does help in understanding the transitions in behavior as a 
function of T [26,29]. The inelastic ball motion, illustrated in fig. 01 models the center of 
mass motion of a granular layer at small T, where the layer remains fairly compact and 
is highly dissipative as a consequence of multiple collisions. We will refer to the inelastic 
ball as "the layer" , meaning the motion of the center of mass of the granular layer. 

For r > 1, the layer leaves the plate at the point in each cycle when the plate 
acceleration exceeds — g. The layer continues in free flight until it later collides with 
the plate. In the regime with squares or stripes oscillating at / /2, the layer leaves and 
hits the plate every cycle, as fig. Eta) illustrates. At T sa 4.0, there is a bifurcation in 
the dynamics, illustrated by fig. Efb): successive trajectories now have different flight 
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Fig. 3. - Phase diagram for granular patterns observed in a vertically oscillated container, 
as a function of the dimensionless acceleration V and dimensionless frequency /* = ,f\/h/g. 
The transitions from a flat layer to squares are hysteretic: solid lines denote the transition for 
increasing T while dotted lines denote decreasing T. (Bronze spheres, a = 0.165 mm; layer 
depth, 5. Oct; container diameter, 770ct.) From [29]. 



durations, with a short trajectory initiated by a collision of the layer with the plate 
whose acceleration greater than g, followed by a long trajectory initiated by a collision of 
the layer with the plate whose acceleration is less than g. This bifurcation corresponds 
to a value of T close to that for the onset of hexagonal spatial patterns [26] , pictured in 
fig-GIc). While square or stripe patterns have the same appearance whether an image is 
obtained at a time t or t + r, hexagonal patterns at these two times are different. The 
two phases are both present and are separated by a phase discontinuity in fig. Hfc) - the 
pattern on the right consists of a hexagonal array of dots, while the hexagonal pattern 
on the left is cellular; the two patterns will be switched one period, r, later. 

Another bifurcation in the dynamics occurs for T > 4.5, where the layer flight dura- 
tion exceeds r, and the layer hits the container every other cycle. Now the velocity of the 
layer relative to the plate at the instant of collision goes to zero (at about T = 4.6), and 
the layer makes a soft landing; not enough momentum is transferred from the vertical 
to horizontal direction to form patterns. When T is increased above about 5.4, there is 
again a transition from a flat layer to a pattern, stripes at low frequencies and squares 
at high frequencies. However, now the pattern period is //4 instead of f/2. Further 
increase in T leads to another bifurcation from trajectories with a single period to sue- 
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Fig. 4. - Trajectory of a completely inelastic ball on an oscillating plate. This is a model for the 
motion of the center of mass of a granular layer. The sinusoidal curve is the trajectory of the 
plate. The ball leaves the plate when the acceleration of the plate becomes -g, that is, when the 
dot-dashed line intersects the trajectory of the ball. If the ball collides with the plate above the 
dot-dashed line, as in (b) and (d), it leaves the plate immediately. From [29]. 



cessive trajectories with different flight durations (fig. Efd)). Again this bifurcation in 
the dynamics of a ball corresponds to the bifurcation of the patterns from squares or 
stripes to hexagons. 

At larger V the layer becomes sufficiently dilated so that the inelastic single ball 
model no longer provides a useful description of the dynamics. But even the //3 and 
//6 regimes of the single ball model (8 < T < 11) have been observed transiently in 
molecular dynamics simulations and laboratory experiments [29] . 

3. Molecular Dynamics simulations 

Simulations of granular media using rigid particles, soft particles, and Monte Carlo 
methods have been conducted by many researchers since the 1980s [7,15,31]. We con- 
sider a granular layer modeled as a collection of hard spheres that interact only through 
instantaneous binary collisions. Between successive collisions the particles move only 
under the influence of gravity. This is known as an Event Driven (ED) type of Molecular 
Dynamics (MD) simulation [32]. Linear and angular momentum are conserved in colli- 
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sions, while energy is dissipated through collisions and surface friction. The interaction 
between particles is described by the coefficient of restitution e, the coefficient of sliding 
friction /i, and the rotation coefficient of restitution f3. (The same values of e, /x and (3 are 
used for ball- wall collisions as for ball-ball collisions.) This collision model was developed 
by Walton [33] and used in hard sphere MD simulations conducted by Bizon [25] and 
Moon [29], as we shall describe. 

The restitution coefficient is often taken to be a constant in MD simulations, but a 
constant value of e can lead to "inelastic collapse" , where particles undergo an infinite 
number of collisions in a finite amount of time [34-36]. However, physically e must 
approach unity as the relative normal velocity v n approaches zero (the elastic limit). 
The simulations presented here assume a form for e that goes to unity in the v n — ► 
limit: e(v n ) = 1 — Bv n ^ for v n < vq, where vq is a crossover velocity, and e(v n ) = eo = 
constant for v n > vq. The approach of e toward unity for low v n avoids inelastic collapse, 
while the constant value eo at high v n is computationally efficient. The MD results are 
insensitive to the form assumed for e(v n ) for v n < vq, as long as e(v n ) increases to unity 
as v n approaches zero [25]. 

In a collision the tangential impulse is given by \i times the normal impulse, with a 
cutoff corresponding to the crossover from a sliding contact to a rolling contact. The 
crossover ratio of the relative surface velocity after collision to that before the collision 
is given by /?, which we fix at -0.35, as in [33]. The parameters eo and fx are set to 
0.7 and 0.5, respectively, values obtained by fitting MD simulation results to laboratory 
observations for square patterns in layers of lead particles at one set of conditions [25] . 

Results from simulation and experiment are compared in fig. for seven values of 
(/*, T). At every point in the phase diagram in fig. 01 the results from the MD simulation 
agree with experiment, not only in the form of the pattern but also in the wavelength of 
the pattern. The dispersion relation relating the wavelength A to /* for layers of varying 
depths collapse onto a single curve, X/h = 1+1. 1/* -1 ' 32 °' 03 , when the container velocity 
exceeds a critical value, v gm ~ 2>^/ag, where v gm corresponds to a transition in the grain 
mobility (gm): for v > v gm , there is a hydrodynamic-like horizontal sloshing motion of 
the layer, while for v < v gm , the grains are essentially immoblile and the stripe pattern 
apparently arises from a bending of the granular layer [24] . 

Most MD simulations and theoretical analyses of granular media consider frictionless 
spherical particles [37-39]. However, MD simulations comparing the behavior of fric- 
tional and frictionless particles indicate that the effect of friction cannot be mimicked by 
increasing the dissipation (decreasing e) [23]; thus friction is not merely an additional 
mechanism of dissipation. Even a small amount of friction increases the overall dissipa- 
tion significantly, not because the frictional dissipation is significant in each collision, but 
because the friction reduces the grain mobility and increases the overall collision rate. 
MD simulations with frictional particles yield square and hexagonal patterns like those 
observed in experiments (fig. [SJ , while simulations without friction do not yield square or 
hexagonal patterns, even if the restitution coefficient is decreased to compensate for the 
absence of friction [23]. Simulations of frictionless particles do yield stripe patterns, but 
the critical T for the onset of patterns is smaller (T c s» 1.9) than for frictional particles 
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Fig. 5. - Standing wave granular patterns in laboratory experiments and molecular dynamics 



simulations for the same number of particles, 60,000, which fill a 100cr x 100(7 container to a 
depth of 5.4 layers: (a) squares, (b) stripes, (c) and (d) alternating phases of hexagons, (e) flat 
layers, (f) squares, (g) stripes, and (h) hexagons. Patterns (a)-(e) oscillate at f/2, (f)-(h) at 
//4. The experiment used lead particles with a — 0.55 mm. From [25]. 



(r c ss 2.5), and the stripes formed by frictionless particles are less robust than those 
formed by particles with friction [23]. 

4. Localized structures and lattice dynamics 

Localized stable standing wave structures dubbed "oscillons" can occur in oscillating 
granular layers when T is decreased slightly below the value corresponding to the onset of 
squares with increasing T [28,40]. Top and side views of oscillons are shown in fig. El An 
oscillon is a small, circularly symmetric excitation that oscillates at f/2; during one cycle 
of the container, it is a peak; on the next cycle it is a crater. Unlike solitons, oscillons 
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Fig. 6. - (a) Snapshot of a container with two oscillons (viewed from above); one is a peak 
(upper left) and the is other a crater (near the center), (b) and (d), Close-ups of an oscillon 
crater viewed from the top and from the side, respectively, (c) and (e) Close-ups of an oscillon 
peak viewed from the top and from the side, respectively. Individual bronze spheres (a — 0.165 
mm) are discernible in (b)-(e). (/ = 25 Hz, T =2.45, layer depth h = 17a.) From [28]. 



are stationary (nonpropagating). Oscillons form with equal probability at all locations 
in the container, and they live indefinitely. If the container acceleration were increased 
slowly from rest to a value just below the onset of squares, no oscillons would appear, 
but an oscillon can be formed by a finite amplitude perturbation (a puff of air or a poke 
with a rod). When oscillons are obtained by decreasing T from the regime with square 
patterns, the number of oscillons that form is not fixed; as many as fifty oscillons were 
observed in the experiment in fig. but no oscillons occurred if T was quasi-statically 
decreased [41]. 

Oscillons of like phase show a repulsive interaction that has a range not much larger 
than the diameter of an oscillon, while oscillons of opposite phase that are closer than 
about 1.4 oscillon diameters attract and form a stable dipole structure, as shown in 
fig. Eta) [28,40], and more complex structures like the tetramer pictured in fig. 0b) and 
the polymer chain in pictured in fig. EI )- 

Can oscillons be considered as building blocks ("atoms") of the square lattice that 
forms with an increase in F? This view is suggested by the observation of the formation 
of a square lattice as T is slowly increased (see fig. Eld) ) : an oscillon seeds a square lattice 
by spawning oscillons, which adjust to form a square array. This observation suggests 
that the granular lattice could be modeled as a system of coupled oscillon atoms, each of 
which is comprised of hundreds of particles that are colliding hundreds of times during 
each oscillation cycle. Thus the lattice approach is much simpler than the full description 
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Fig. 7. - (a) Dimer formed by two bound oscillons of opposite phase; one period of oscillation 
of the container later the white peak will have become a crater and the black will have become 
a peak, (b) Tetramer formed of four oscillons. (c) Polymer chain of five oscillons. (d) A square 
lattice grows by nucleating oscillons. Individual bronze spheres (a — 0.165 mm) are discernible 
in (a)-(c). From [28] 



of all the particle motions and simpler than a continuum fluid description of the granular 
medium. As we will now describe, the lattice picture is supported by an analysis of the 
dynamics of the lattice. 

Close examination of the center of mass of a peak in a square lattice reveals harmonic 
motion about the equilibrium position of the peak for a wide range of T and / [27]; such 
a lattice oscillation is illustrated in fig. |SJ One test of the conjecture that the lattice 
of peaks can be modeled by a lattice of balls connected with springs is to compare the 
dispersion relation for the two lattices. A time sequence of images of the granular pattern 
was Fourier-transformed in space and then in time to obtain the frequencies of oscillations 
of lattice modes with different wave vectors. The frequency /l of the lattice modes as a 
function of k (the magnitude of the wave vector) was found to be well described by the 
dispersion relation for balls connected by springs, Jl = fBz\sin(ka/ (2y/2)\, where fsz is 
the frequency at the edge of the Brillouin zone and a the lattice spacing [27] . The lattice 
oscillation frequency _/x is typically an order of magnitude smaller than / but depends 
on the plate acceleration T and frequency /. 

In a crystalline solid, defects form when the amplitude of oscillation of atoms about 
their equilibrium position in the lattice becomes large. To investigate the possible forma- 
tion of defects in a granular lattice, the plate frequency / was modulated at the lattice 
frequency f^. For sufficiently large modulation amplitude, defects formed, breaking the 
long range order of the square lattice. It was found that the amplitude of the lattice 
oscillations could be increased further by adding a lubricant (graphite powder) to reduce 
the friction between particles. The result was that the granular lattice melted [27]: the 
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Fig. 8. - Left: close up snapshot of a granular lattice (F" = 2.9, / = 25 Hz, h = 4<r). Right: time 
evolution of the peaks in the boxes A and B in the left-hand image. The peaks oscillate out of 
phase with a frequency about twenty times smaller than / — 1/r. (from [27]). 
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spatial Fourier transform became a circular ring (about k — 0) rather than sharp peaks, 
as fig. El illustrates. 

Lattice oscillations and defect formation have also been studied in MD simulations, 
where friction can be easily varied or set to zero. First a square lattice was simulated with 
(i = 0.5, as described previously and illustrated in the first panel of fig- Efb). Then /i was 
set to zero, and defects were observed to form quickly as the lattice oscillation amplitude 
increased. Finally the lattice melted, as illustrated by the last panel of fig. Efb) . 

Studies of melting in two-dimensional solids have shown that melting occurs when the 
Lindemann ratio, 7 = (\u m — u n \ 2 )/a 2 , exceeds 0.1 [42,43]. Here u corresponds to the 
displacements of atoms from equilibrium lattice sites, a is the lattice constant, and the 
average is taken over all nearest neighbors m and n. The Lindemann ratio was computed 
for the granular lattice in the MD simulation, and it was found that when the friction 
coefficient fi was decreased, 7 increased. Further, when 7 reached the value 0.1 (which 
happened for [i — 0.1 for the conditions of the simulation), the granular lattice melted, 
in accord with the result for crystalline solids [27]. 



5. Continuum Description 

Section 3 introduced Molecular Dynamics simulations as a useful tool in describ- 
ing granular flows. This technique models the system on a microscopic level, evolving 
individual particle trajectories using Newton's laws and computing the effects of each 
collision. Averaging over many collisions and particle trajectories gives the macroscopic 
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Fig. 9. - Defect creation and melting of a square granular pattern in a vertically oscillating 
granular layer: (a) experiment and (b) molecular dynamics simulation. The insets show Fourier 
transforms of the spatial patterns. In the experiment the plate oscillation frequency (/ = 32 
Hz) was modulated at the natural oscillation frequency of the lattice (2 Hz), and at t = r 
graphite powder was added to the layer of bronze spheres. By t = 56r defects had formed, and 
by t = 175r the lattice had melted. In the MD simulation the friction coefficient was reduced 
from n = 0.5 to zero at t = r; by t — 22r defects had formed, and by t = lOOr the lattice had 
melted, (r = 2.9.) From [27]. 



behavior of the flow. A complementary method for understanding granular flows is to 
model the macroscopic motion directly by a continuum field theory that describes the 
bulk motion of the flow in terms of the density, velocity and temperature fields. Unlike 
MD simulations, the continuum approach is not limited by particle number. A personal 
computer currently contains enough memory for useful MD simulations of laboratory 
experiments. However, industrial processes contain billions of particles, far outside the 
abilities of MD simulations. Another reason that a continuum approach is attractive is 
that it could exploit tools such as stability analysis, amplitude equations, and perturba- 
tion theory, which have been developed through more than a century of research on the 
the Navier-Stokes equations and other partial differential equations. 

Granular flows present many difficulties in developing a continuum theory [44,45]. 
Continuum theory requires a separation of length and time scales: variations over space 
should be small and occur over long distances, so that the behavior of local collections of 
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Oscillated Granular Layer Oscillated Liquid Layer Ray leigrhBernard Conwction 



Fig. 10. - Forced granular materials produce qualitatively similar patterns as forced fluids: (a) 
stripe pattern formed by a vertically oscillated granular layer [26] , (b) stripe pattern formed by 
a vertically oscillated layer of water [52], (c) stripe pattern formed in thermal convection of a 
fluid (COa) [55]. 



individual particles can be averaged and replaced with small fluid elements. Changes in 
time for the flow should occur for times long compared to the mean time between particle 
collision so that particles moving between fluid elements do not affect the average values 
in a fluid element. Unfortunately, inelastic collisions between particles create an inherent 
lack of scale separation [11,45]. Sufficient separation of scales may only be present for 
granular flows in the specific circumstances of low density and low dissipation [11,45,46]. 

The derivation of the continuum equations from kinetic theory makes assumptions 
about the underlying statistics of granular flows, assumptions which have not been ver- 
ified by MD simulations. For instance, the velocity distribution function is assumed 
to have a steady state functional form that is nearly Gaussian. Since granular flows 
are dissipative, a steady state distribution function can only be achieved in the pres- 
ence of forcing. Granular experiments have yielded velocity distributions that depend 
on the forcing characteristics and experimental geometry [47-50]. Also, most deriva- 
tions of continuum equations assume Boltzmann's molecular chaos (particle velocities 
before collisions are uncorrelated) , but strong velocity correlations have been found in 
MD simulations [51]. 

Despite the reservations regarding a continuum approach in granular media, obser- 
vations of granular media have revealed many phenomena similar to those observed in 
continuum systems. For example, the stripe patterns shown in fig. llOf a') look like those in 
vertically oscillated liquid layers [52] ffig. lluf b)). chemical reaction-diffusion systems [53], 
Rayleigh-Benard convection in fluids [54] (fig. HUT c). and liquid crystals [56]. 

Not only are the patterns similar for granular and continuum systems, but also some 
the same pattern instabilities have been observed. For example, when the wavenumber 
of parallel convection rolls (stripes) in a Rayleigh-Benard convection becomes small, 
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Cross Roll Skew- Varicose Spiral Defect 
Instability Instability Chaos 




Fig. 11. - Instabilities of patterns found in oscillating granular layers and Rayleigh-Benard 
convection in a fluid. Cross roll instability in stripes: (A) Vibrated granular layer [57] and (a) 
Rayleigh-Benard convection [59]. Skew varicose instability in stripes: (B) Granular layer [57] 
and (b) Rayleigh-Benard convection [60]. Spiral defect chaos in: (C) vibrated granular layer [61] 
and (c) Rayleigh-Benard convection [60]. 



an instability leads to the formation of cross rolls with a larger wavenumber that are 
perpendicular to the original rolls [57,58]; the same instability has been observed for 
stripes in oscillated granular layers, as fig. ^2 illustrates. The cross rolls invade the 
region of small wavenumber stripes such that, after sufficient time, the region contains a 
pattern of straight stripes perpendicular to the original pattern and with a larger wave 
number. 

Granular stripe patterns also exhibit a skew varicose instability like that in convection 
roll patterns ffig. lll|) . When the local wavenumber becomes too large, an initially straight 
pattern of stripes will develop a distortion which evolves into a dislocation defect. The 
defect propagates away; destroying one of the stripes and decreasing the local wave num- 
ber of the pattern. The stability of the stripe pattern in fluid convection is well described 
by amplitude equations derived from the Navier-Stokes equations for fluids [58]. That 
the granular pattern shows the same behavior strongly suggests a continuum description 
for the vibrated system is applicable. 

Aspects of the phase diagram for granular patterns (fig. 13) have been reproduced by 
amplitude equation models. For example, a phenomenological continuum model requiring 
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that the mass of the layer is conserved locally, produces stripe, square, and oscillon-like 
patterns similar to those found in experiment [62,63]. A continuum, shallow water like 
model of the granular layer captures the patterns and yields a dispersion relation which 
agrees with experiment [64]. The success of these and other models [65-68] provides 
further motivation for considering continuum equations derived for a granular gas. 

Additional evidence for the applicability of continuum theory to granular media is pro- 
vided by a recent study of noise in vertically oscillating granular layers. In the Rayleigh- 
Benard system below the onset of convection, thermal noise has been found to drive 
noisy transient disordered waves with a characteristic length scale. The intensity and 
coherence of these modes increases as the transition from conduction to convection is 
approached [69]. This behavior is well-described by the fluctuating hydrodynamic the- 
ory of Swift and Hohenberg [70]. Remarkably, the same noisy incoherent modes are 
observed just below the transition from a flat vertically oscillating granular layer to a 
square pattern [71]. The Swift- Hohenberg continuum theory describes the observations 
for the granular system very well, even though the noise is not thermal noise, which is 
many orders of magnitude too small; apparently the noise arises from the fluctuations 
due to the small number of particles [71]. 

Fired by the promise of quantitative predictive power and encouraged by the quali- 
tative similarity of granular flows to fluid flows, researchers have proposed various con- 
tinuum descriptions for rapid granular flows [72-78]. This section focuses on one such 
description [75] and compares results from it to MD simulations and a granular flow 
experiment. 

Jenkins and Richmann derived a set of inelastic continuum equations in a manner 
similar to the derivation of the Navier-Stokes equations [75] . This approach begins with 
the single particle distribution function /W (r,£), which gives the probability of finding 
a particle at a position r with a velocity v at a given time t. Integrating /W over all 
possible velocities gives the local number density, n(r, t). The ensemble averaged value 
of any particle property ip is determined by 



(5.1) <^>=-U{v)f^{v,v,t)dv, 

nj 

The Boltzmann equation describes how /W changes in time. Particles can move in 
and out of volume elements due to streaming motion; particle velocities can change in 
response to external forces F; or particles can be scattered out of elements by collisions. 
The time rate of change for /W is given by 



(5.2) %— + v • V r /« + F • V v /« = Q(f^), 

at 

where Q(f^) is the collision operator. Collisions are considered to be binary, frictionless, 
and inelastic with a constant coefficient of restitution eo [75] . Integrating eq. 115.211 yields 
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the balance law for the number density, 



(5.3) 



dn 
~dt 



+ V • (nu) = 0, 



where u(r, t) = (1/n)/ v/W (r, v, t)dv is the local average velocity. Multiplying by the 
velocity and then integrating gives the balance law for momentum, 



(5.4) 



du 

n[ — + u • Vu ) = V • P - ngz, 



Finally, multiplying by v 2 and integrating gives the balance law for the energy, where 
the granular temperature T is proportional to the average kinetic energy of the random 
motion of particles, 



(5.5) 



T= 1/3 (< v 2 > - < v > 2 ) 



(5.6) 



Hf— ) 



-V-q + P :E-7. 



The granular temperature T is many orders of magnitude greater than the Boltzmann 
temperature: thermal fluctuations are negligible (mga » fcsTs). 

A series of approximations is required in order to derive the form of the pressure 
tensor P, the velocity gradient tensor E, and the heat flux q. One assumes that 
is nearly Gaussian, that spatial derivatives of n, u, and T are small, and that (1 — eo) 
is small. With these assumptions, the components of the velocity gradient tensor E are 
given by: — | {djU% + diUj). The components of the stress tensor P are given by the 
constitutive relation: 



(5.7) 



P — 



-P + ( A - ^) E kk 



Sij + 2/j.Eij, 



and the heat flux is given by Fourier's law: 

(5.8) q = -kVT. 

The transport coefficients are fully determined and are the same as for a dense gas of 
hard spheres. The bulk viscosity is given by 



(5.9) 



A 



3V^ 
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the shear viscosity by 



(5.10) 



/' 



6 



and the thermal conductivity by 
15v^ 
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(5.11) t 

where 

(5.12) G{u) = vg Q {v), 

and the radial distribution function at contact, is [79]: 



(5.13) 



50 (f) 



1 - 



where v is the volume fraction of the flow and v„ 



0.65 is the 3-dimensional random 



close-packed volume fraction. 

The only difference between these equations and those for an elastic gas is 7 in 
eq. Q5.6I1 . which accounts for the temperature loss due to inelastic collisions: 



12 nT 3 / 2 
(5.14) 7= _^(i_ e 2)!^ G(2/) . 

The system is closed by an equation of state, proposed by Goldshtein et al. in [79], 



(5.15) 



p = nT [1 + 2(1 + e<,)G(i>)] 



Direct experimental verification of the inelastic continuum equations has been slow 
in coming due to the complexity of solving the equations and also due to difficulties in 
finding an appropriate experimental system [7]. The presence of strong gradients in 
granular materials [11,45] adds additional difficulty to solving continuum equations. For 
instance, simulations of the vertically vibrated layer find that the temperature varies by 
three orders of magnitude throughout the cycle [80]. Thus, unlike most Navier-Stokcs 
simulations, the transport coefficients (A, fj, and k) cannot be treated as constants, but 
must be recomputed at each grid point at every time step. Additionally, a complete set of 
boundary conditions for granular flows is still not established and this remains an active 
area of research [81-86]. Without the correct boundary conditions, numerical solutions 
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can be unstable and are not guaranteed to converge to a correct solution in the bulk. 
For a good discussion on the difficulties in determining the correct boundary conditions, 
see Goldhirsch's review paper [11]. 

In a 1990 review paper [7], Campbell made a resounding call for granular flow exper- 
iments to make quantitative tests of the inelastic continuum approach. The application 
of new technologies such as particle tracking in two and three dimensions is now making 
these measurements feasible, but granular flow experiments still present technical chal- 
lenges. Plugs develop in pipe flow [87,88], wall effects dominate in quasi- two-dimensional 
experiments [89] , and detailed bulk flow measurements are difficult to make in fully three 
dimensional experiments [90]. 

The distinguishing feature of granular flows is that inelastic collisions dissipate energy. 
Without an external source of energy, the granular temperature decays to zero with all 
particles coming to rest. Experimentally, energy can only be put into the flow through the 
boundaries. Shock waves serve as a mechanism to deliver energy from the boundary to 
the bulk of the flow. Studying the balance between the transfer of energy by shock waves 
and the energy dissipation through inelastic collisions is important in understanding 
granular flows [12]. In the next section we present two studies of shock waves as a test 
for the inelastic continuum equations. 

6. Shock Waves in Granular Materials 

The sound speed (the speed a pressure wave travels) in a granular medium is typically 
much smaller than the streaming velocity. Hence shocks are common in granular media. 
For example, imagine pouring sand out of a bucket. Gravity accelerates the flow down- 
ward, creating an average velocity U that reaches 100 cm/s after the sand has fallen by 
only 5 cm. In contrast, the sound speed c in the granular gas becomes small, typically 
10 cm/s, as multiple particle collisions cool the gas, reducing the random velocities of 
the particles. The simple act of turning over a bucket full of sand can easily generate 
a supersonic flow with Mach numbcr=-^=10. (A flow with Mach number greater than 
unity is supersonic.) 

The sound speed in a granular gas can be determined from thermodynamic relations, 



where c is the speed of a sound, p is the local density, S is the entropy, c p is the specific 
heat at constant pressure, and c v is the specific heat at constant volume. For a dense 
inelastic gas, c is given by [91]: 



(6.1) 




(6.2) 
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where x = 1 + 2(1 + e)G{v). 

A shock forms when a supersonic flow encounters an obstacle. The shock separates 
two regions of the flow, an undisturbed region that is unaware of the obstacle, and a 
compressed region which has adjusted to fit the boundary conditions at the obstacle. The 
compressed region has a higher temperature and smaller velocity than the undisturbed 
region. In an ideal fluid with no viscosity, heat conduction, or dissipation, a shock is a 
zero-width surface of discontinuity. In a non-ideal fluid the shock has a finite width on 
the order of a particle's mean free path in the fluid [92]. 

When a fluid with velocity U > c impinges perpendicularly onto an obstacle, a normal 
shock forms and propagates in the —U direction. Section 6.2 will discuss the formation 
and propagation of normal shocks in a vertically vibrated granular layer where the flow 
fields and thus the shocks are highly time dependent [80]. If, instead, the fluid velocity 
and the obstacle are not perpendicular, an oblique shock forms and propagates into the 
flow at an angle and with a speed determined by the local flow values. Section 6.1 
describes an oblique shock formed in a steady state laboratory flow past a wedge [89] . 

6T. Steady state flow past an obstacle. - We now describe an experimental study of 
shocks in a time-independent flow, and then the experimental results will be compared 
to two simulations: an event-driven molecular dynamics simulation, similar to those 
discussed in Section 3, and a two-dimensional, finite-element solution of the inelastic 
continuum equations presented in Section 5. 

In the experiment, stainless steel spheres (particle diameter a = 1.2 mm) fell under 
gravity past a wedge sandwiched between two glass plates separated by 1.6cr. The par- 
ticles were initially distributed uniformly on a conveyor belt. As the conveyor turned, 
particles fell off into a hopper that guided the particles into the cell formed of the closely 
spaced plates; the wedge was located a distance of 42a below the top of the cell. The 
positions and velocities of the particles were determined from high speed images of the 
falling particles, and data from many thousands of particles were averaged to obtain 
the time-independent velocity, volume fraction, and temperature fields. The average 
free stream speed of sound determined from the measurements for flow incident on the 
wedge was 0.09 m/s. The flow entered the top of the cell with a Mach number of 7 and 
accelerated under gravity to a Mach number of 12 at the tip of the wedge. 

The horizontal velocity field measured in the experiment is shown in fig. I12f a~) . A 
shock separates the undisturbed region, where the horizontal velocity is nearly zero, from 
the compressed region, whose stream lines follow the flow around the obstacle. Because 
of gravity and inelasticity, the shock does not extend out at a constant angle but curves 
towards the wedge. 

At the bottom of the wedge the compressed gas expands in a fan-like structure as 
the volume available to the flow increases (fig. I13|) . In an expansion fan the density and 
temperature decrease and the Mach number increases. The expansion fan is a smooth 
transition radiating from the bottom corner of the wedge. 

The flow was computed numerically in a three-dimensional MD simulation (fig-Elk)) 
and in a two-dimensional finite difference simulation of the inelastic continuum equations 
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Fig. 12. - Horizontal component of the velocity field of a granular flow incident downward on a 
wedge, determined by three methods: (a) experiment, (b) MD simulations, and (c) integration of 
inelastic continuum equations. Each picture shows a region 130<r by 104<j. The solid lines with 
arrows denote streamlines. Quantitative comparisons along the dashed line in (a) are shown in 
figs. E31 and [TE] (From [89]). 



(fig. H2f c)). The two simulations yield results for the horizontal component of velocity in 
qualitative accord with experiment: a shock forms at the tip of the obstacle, and behind 
the shock the flow is compressed, has a higher temperature, and lower mean velocity. 
Quantitative comparisons among the methods are plotted for values of the fields along 
the dashed line shown in fig. I12f a). 

Three parameters were adjusted in the MD simulation to achieve the agreement with 
the experiment shown in figs. El The same coefficient of restitution eo = 0.97 and 
friction coefficient (i — 0.15 were used to model interparticle and particle-wall collisions. 




Fig. 13. - The horizontal velocity field measured for the expansion fan that formed when the 
supersonic granular flow reached the bottom of the wedge. The solid lines indicate selected 
streamlines. The total height of the region shown is 55a". The white region below the wedge 
had too few particles for the velocity to be determined. (From [89]) 
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Fig. 14. - Shock profiles for granular flow past a wedge measured in an experiment (circles) are 
compared with results from molecular dynamics (solid lines): (a) volume fraction, (b) horizontal 
component of the velocity, and (c) temperature. The profiles are taken along the dashed line in 
fig. El (From [89]) 



The initial conditions of the experiment were modeled by placing particles into the top 
of the cell at a constant rate. Incoming particles were placed randomly at the top of 
the cell with a mean downward velocity measured from the experiment, and fluctuations 
were chosen from a Gaussian distribution determined by the measured temperature. An 
additional parameter a, defined as the ratio of temperature perpendicular to the wall to 
that parallel to the wall, was set to 0.8. These parameters, which were not measured in 
the experiment, were adjusted to provide agreement in the full flow fields, including the 
free-stream velocity. 

Results from the MD simulation are compared with experiment in fig. 1141 for the 
volume fraction, horizontal velocity component, and temperature. The agreement is 
quite good with a root mean square difference between experiment and simulation of less 
than 2% for the volume fraction and velocity fields and 10% for the temperature field. 

The simple geometry and steady state behavior of the experiment provided a good 
system for testing the inelastic continuum equations. However, a full three-dimensional 
simulation of the experiment was found to be time prohibitive. Instead, the equations 
were solved on a two-dimensional grid; consequently the simulation could not capture the 
interaction of particles with the confining glass side walls. Frictional collisions with the 
side walls strongly affected the flow in the experiment and in the fully three-dimensional 
MD simulations. The average downward acceleration of a single particle falling between 
the two glass plates in the experiment was 8.9 m/s 2 , while the same particle falling 
outside the cell accelerated with the expected 9.8 m/s 2 . 

The continuum equations were numerically solved by a second-order accurate, finite 
difference method. The only fit parameter in the equations was the coefficient of restitu- 
tion, which was set to the same value of eo(0.97) used in the MD simulation. Boundary 
conditions at the inlet were determined by the experiment and at the outlet were free. 
Slip velocity boundary conditions were used along the wedge boundary. The heat flux at 
the wedge was taken to be proportional to the local v and T 3 / 2 [85]. Euler time stepping 
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Fig. 15. - Comparison of shock profiles for granular flow past a wedge obtained from molecular 
dynamics (solid lines) and inelastic continuum equations (dotted line), assuming no friction, 
(a) Volume fraction, (b) horizontal velocity profile, and (c) temperature along the dashed line 
in fig. Il^f a). Experimental measurements (open circles) show similar qualitative behavior but 
disagree quantitatively. The difference between the simulations and the experiment is due to 
wall friction. (From [89]) 



was used to increment the simulation until the flow reached a steady state where the 
horizontally averaged mass flux was constant to 0.01 

Experiment and continuum simulation showed similar behavior, but the shape of the 
curves differed and the magnitudes of the fields disagreed by as much as a factor of 
two (fig. I15|) . This disagreement was attributed to the factional drag of the confining 
side walls in the experiment. A three-dimensional simulation of the inelastic continuum 
equations with viscous boundary conditions along the side walls should agree better with 
the experiment [12]. 

Molecular dynamics simulations were done with wall drag neglected for comparison 
with the two-dimensional simulation of the continuum equations. The two simulations 
agreed remarkably well in all regions of the flow except within 5a of the wedge tip. 
Near the tip of the wedge, the two simulations disagreed due to different boundary 
conditions. The agreement between the two simulations in the bulk of the flow confirms 
the applicability of the continuum description for granular flows. The disagreement with 
the experiment emphasizes once again the importance of including friction in a continuum 
description, both in the equations and in the boundary conditions. 

6'2. Shock waves in a vertically oscillating layer. - Section 2 described the patterns 
that form when a layer of granular materials is vertically oscillated. Shocks play an 
important role in this system [80]. Each time the layer collides with the plate a shock 
forms and propagates through the layer, transmitting energy upward through the layer 
to the surface. The striking similarity of the granular patterns to their counterparts 
in continuum systems strongly suggests that a continuum description of granular flows 
could prove useful. If this is true, the continuum description must also capture the shock 
dynamics in the bulk of the layer. 

Here we compare inelastic, frictionless, fully three-dimensional MD and continuum 



Pattern formation and shocks in granular gases 



25 







ft=0.lS 










— I 1 -, 1 ii ii 1 region 


ptatR 




ft=0.22 








TV 


plate 



10- 







I 




I'- 


plate 


ll 


ft=0.40 



0.2 0.4 
v, (5G/ga)T 



0.6 



Fig. 16. - Dimensionless temperature T/ga and volume traction v as functions of the dimen- 
sionless height z/a at four times ft in the oscillation cycle. For each time, a snap shot from 
the MD simulation is shown in the left column, with individual particles color coded according 
to temperature: high T in red, low T in blue, and the bottom plate of the container shaded 
solid gray. The right column shows horizontally averaged v (blue) and T/ga (red) for the same 
four times. The plate is shown as a horizontal black solid line, results from MD simulation are 
shown as dots, and continuum results are solid lines (From [80]). 
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simulations. In order to focus on the formation and propagation of the shock wave, 
pattern formation is intentionally suppressed by considering a container smaller than 
one wavelength of the pattern in either horizontal direction. We focus on one point in 
the phase space shown in fig. [3] T — 3 and /* = 0.42 for a layer with depth h — 9a. 
For this set of parameter values, a vibrated layer in a larger cell would have a //2 stripe 
pattern. 

In both simulations periodic boundary conditions were used in the two horizontal 
directions and impermeable boundary conditions, u z = 0, were applied at the plate. 
The additional boundary conditions required for the continuum simulations were taken 
from the MD simulation. In the MD simulation, the vertical derivatives at the plate 
were negligible throughout most of the cycle. For simplicity, the continuum simulation 
assumed du x /dz = 0, du y /dz = 0, and dT/dz = at all times in the cycle. 

The evolution of the shock wave throughout a plate cycle is shown in figs.ED an d IT7I 
The dynamics of the cycle occurs in the time interval between ft = and one cycle later, 
ft = l. 

At ft = the container is at its minimum height. The layer, having been thrown off 
the plate in the previous cycle, now falls towards the plate. Inelasticity has dissipated 
most of the energy so that the layer's temperature is nearly zero. The Mach number of 
the layer with respect to the plate is much greater than one. The MD and continuum 
simulation show similar behavior in the v and T fields. 

At ft = 15 the layer begins to collide with the plate. A shock wave forms, separating 
the region near the plate where v and T increases from the undisturbed region still falling 
towards the plate. 

At ft — 0.22 the shock wave is moving through the layer. The compressed region 
continues to grow. Collisions between particles in this high density region cause the layer 
to cool behind the shock, creating a lower temperature near the plate. 

At ft = 0.40 the shock has propagated through the layer and into the very dilute 
region above the layer. At this time, the plate is approaching its maximum height and 
the layer begins to leave the plate as the downward plate acceleration exceeds g. The 
layer continues to cool behind the shock, setting the stage for the next oscillation. 

The MD and continuum simulations show good agreement throughout the cycle, de- 
spite the presence of large spatial gradients and a strong time dependence. In the dilute 
regimes above and below the layer, numerical solutions of the inelastic continuum equa- 
tions are unstable unless artificial dissipation is added [80], following the example from 
numerical solutions of Knudsen gases. [93]. The effect of the extra dissipation is most 
pronounced in the falling layer and accounts for the disagreement between MD and con- 
tinuum at the top of the cell. 

Inelastic collisions between particles are the distinguishing characteristics of a granular 
gas, but few studies have examined how granular flow properties depend on the restitution 
coefficient. Simulations for the oscillating layer were modified to study the effect of 
varying eo on the propagation of the shock. The initial conditions for this numerical 
experiment were taken for a layer with eo = 0.99, using the same parameters as in the 
above discussion. At ft = 0.33, when the center of mass layer was near its maximum 
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Fig. 17. - Location of the shock (solid line for continuum, squares for MD) and the center of 
mass of the layer (dashed line for continuum, circles for MD) as a function of time ft during one 
cycle of the plate (thick solid line) for particles with eo = 0.90. The plot is shaded according to 
the volume fraction from the continuum simulation, so that high volume fraction is dark and low 
volume fraction is light. The "top" and the "bottom" of the layer from MD (when the volume 
fraction drops to less than 4% of that for random close packed particles) are shown as +'s and 
X's respectively. The material below the shock is compressed as compared to the region above 
the shock, as can be seen from the shading. From( [80]). 



height above the plate, the coefficient of restitution was suddenly changed, which changed 
the subsequent evolution of the shock. 

As before, when the layer hits the plate it compresses and forms a shock that prop- 
agates through the layer. The smaller the value of eo, the faster the layer cools and 
compacts; for small eo, the layer remains very compact throughout the cycle and leaves 
the plate almost as a solid body. For higher values of eo, the layer dilates quickly after 
each collision with the plate. The maximum height of the center of mass in a cycle 
increases with increasing eo- 

The speed of the shock (eq. (JBJ) depends on both the temperature of the flow and on 
its density. Since the density and temperature of the flow change throughout the cycle, 
so does the shock speed as the shock propagates through the layer. The behavior of the 
average speed of the shock as a function of inelasticity is shown in fig. ^| For small 
values of eo the shock speed asymptotes to a fixed value of The shock speed 

monotonically increases with increasing eo- The special case of elastic particles appears 
to match with the limit of e — > 1, suggesting that there is no qualitative difference in 
shock propagation through elastic and inelastic gases. 

For both problems we have considered - granular flow past a wedge and a vertically 
oscillating granular layer - numerical solutions of the inelastic continuum equations of 
Jenkins and Richman agree well with MD simulations for frictionless particles. The con- 
tinuum equations were derived for a weakly dissipative, low density, frictionless granular 
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Fig. 18. - Average dimensionless shock speed, v s h oc k/ \f§®, in the reference frame of the plate. 
Vshock is calculated as the average speed of the shock from when the shock is formed until it 
leaves the layer. (From [80]) 



gas, assuming small gradients in the flow fields in both space and time. Nevertheless, the 
equations capture the evolution of a shock through a dense, inelastic oscillating layer, 
and qualitatively capture the properties of a shock formed in flow past a wedge. With 
more research on boundary conditions and the incorporation of friction, these continuum 
equations show great promise. 

7. Discussion 

More than one thousand papers have been published on granular materials since de 
Gennes brought the subject to the attention of physicists, and Bak's work (1987) on 
self-organized criticality stimulated interest in sand piles [94]. However, much remains 
to be done to achieve a level of understanding of granular media comparable to that for 
fluids and crystalline solids. Experiments and simulations have investigated a wide range 
of problems including the angle of repose [95] and internal structure of sand piles [96,97], 
shear forces in Couette- Taylor flows [90,98], convection due to temperature gradients [99- 
101] and due to buoyancy [102], flows in a rotating drum [103,104], and chute flows [31, 
105-107]. Much of the research has concerned granular media as a solid where particles 
are in continuous contact, while this chapter has concerned rapid granular flows (the 
"collisional regime" ) where inertial effects are important and force chains do not play a 
major role. We have further limited the considerations to particles that interact only on 
contact. 

Understanding flows of grains that interact only on contact would seem at first to 
involve a straightforward application of Newton's laws. However, the energy loss in 
collisions complicates the application of standard statistical methods. We have focused 
on two systems that hint at the rich variety of phenomena exhibited by granular media in 
the collisional regime: a vertically vibrated granular layer and supersonic granular flow 
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past an obstacle. These two problems were chosen because they are amenable to direct 
comparison of experiment, molecular dynamics simulations, and continuum theory. 

A vertically vibrated granular layer spontaneously forms spatially extended patterns. 
Experiments and MD simulations reveal that the collective motion of grains arises due 
to dissipative collisions between particles, and does not require mediation by an intersti- 
tial gas or side walls. For square granular patterns, an approach intermediate between 
molecular dynamics and continuum models has been found to describe the dynamics of 
the lattice pattern: a collection of particles that form a peak (an oscillon) is like an atom 
in a crystalline lattice. The modes of the granular lattice obey the dispersion relation for 
a two-dimensional lattice, and the granular lattice even forms defects and melts in the 
same way as a two-dimensional crystal of atoms (Section 2). 

The spatial patterns formed by oscillating granular layers exhibit marked similarities 
to those observed in continuum nonequilibrium systems such as convecting fluids and 
oscillating liquids (Section 5). Further, the cross-roll and skew- varicose instabilities ob- 
served in thermal convection in a fluid and interpreted in terms of the hydrodynamic 
equations (more specifically, the Boussinesq equations) have also been observed in oscil- 
lated granular layers. Various amplitude equation models have been found to describe 
granular patterns and their instabilities. Even the subtle effects of noise on the transition 
from conduction to convection in fluids have been found also in oscillating granular layers 
near the onset of the transition from a flat layer to a square pattern. 

The striking similarities of granular patterns to those found in nonequilibrium con- 
tinuum systems and in experiments on granular flows under shear and in rotating drums 
suggest that granular gases may be describable by continuum theory. Inspired by these 
observations, researchers have proposed many continuum descriptions. The descriptions 
differ in the particle properties included in the collision model; for instance, collision 
models can be frictionless [75,108] or can account for friction between particles [109,110]. 
Equations of motion obtained by Goldshtein and Shapiro include, in addition to the terms 
in the Navier-Stokes equation, a term accounting for heat transport by density gradi- 
ents [76]. A Model presented by Bocquet et al. [90] includes corrections to the viscosity 
due to velocity correlations. None of these models has been definitively established. 

We have compared predictions of continuum equations derived by Jenkins and Rich- 
man with experiments on shocks in vibrating layers and flow past an obstacle. For 
both geometries, numerical solutions of the inelastic continuum equations agree well 
with results from MD simulations of smooth (frictionless) inelastic spheres. However, 
comparisons of the continuum equations with experiment and with MD simulations for 
particles with friction have demonstrated the crucial role of friction in granular flows. 
For continuum equations to achieve quantitative predictive power, the effects of friction 
between particles and between particles and boundaries must be included. 

Derivations of granular hydrodynamic equations have thus far assumed weak dissi- 
pation and small particle volume fractions. Future work should extend theory to higher 
densities and larger dissipation. Because inelastic collisions dissipate temperature, gran- 
ular flows frequently coexist with solid-like sand piles. A major challenge is the develop- 
ment of theory that bridges the gap from the collisional regime to the quasi-static regime 
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where particles are always in contact. 

As continuum equations of motion become better established, it will become possible 
to exploit the power of the continuum description. Continuum models for larger and 
denser systems may reveal new phenomena. Continuum simulations are better suited 
to time-dependent flows than MD simulations. Linear and nonlinear stability analyses 
could provide insight into bifurcations, just as a century of stability analyses of the 
Navier-Stokes equations has given insights into diverse fluid flow phenomena. Stability 
analyses have been conducted for simplified continuum models (see, e.g., [Ill]), but thus 
far no stability analysis has been conducted for a realistic set of granular hydrodynamic 
equations. Further, given a set of granular equations like the Navier-Stokes equations, it 
should be possible to derive amplitude equations, which can yield a better understanding 
of instabilities in granular flows. 

The experiments and theory presented in this chapter involve spheres of uniform size, 
while industrial applications usually involve a wide range of particle sizes and shapes. 
Experiments and simulations on flows with two particle sizes show an additional richness 
to granular flow phenomena such as size segregation [23,112-115], nonequipartition of 
energy [116,117], and increased normal stress [117]. Much theoretical and experimental 
work is needed on systems of particles with a range of sizes and shapes. 

Air friction is usually neglected in simulations and in the interpretation of experi- 
ments, whether or not the experiments are conducted in vacuum. In contrast, in gran- 
ular systems in industry, air friction and buoyancy are often important, although the 
air effects can be negligible for large particles (say greater than 1 mm). The interaction 
of the interstitial fluid and sand leads to the development of sand dunes [118] and sand 
ripples [119, 120] and the formation of heaps in vibrated layers [121]. The inclusion of in- 
terstitial flow in continuum theories for granular materials is another challenge for future 
research on granular materials. 

In conclusion, much remains to be done to establish a fundamental understanding 
granular flows. Future studies should seek more examples of granular flows amenable 
to experiment, molecular dynamics simulations, and continuum theory. A concerted at- 
tack from three approaches should lead to a better understanding of outstanding issues 
concerning boundary conditions, the incorporation of particle friction and velocity cor- 
relations into theory, and the extension of theory to higher particle volume fractions and 
higher dissipation. 
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